/*
 * compute the SPHARM coefficients and the associated Mesh
 *
 * author:  Martin Styner
 *
 */

#ifdef _MSC_VER
#pragma warning ( disable : 4786 )
#pragma warning ( disable : 4503 )
#pragma warning ( disable: 4284 )
#endif

#include <time.h> // by feilang
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <string>
#include <string.h>

#include <itkMeshSpatialObject.h>
#include <itkMesh.h>
#include <itkSpatialObjectWriter.h>
#include <itkSpatialObjectReader.h>
#include <itkTriangleCell.h>
#include <itkDefaultDynamicMeshTraits.h>
#include <itkProcessObject.h>

#include "ParametricMeshToSPHARMSpatialObjectFilter.h"
#include "SphericalHarmonicSpatialObject.h"
#include "SphericalHarmonicMeshSource.h"
#include "SphericalHarmonicCoefficientFileWriter.h"
#include "SphericalHarmonicCoefficientFileReader.h"
#include "itkMesh3DProcrustesAlignFilter.h"

#include "vtkPolyDataToitkMesh.h"
#include "itkMeshTovtkPolyData.h"
#include "vtkPolyDataReader.h"
#include "vtkPolyDataWriter.h"
#include "vtkIterativeClosestPointTransform.h"
#include "vtkLandmarkTransform.h"
#include "vtkPointLocator.h"
#include <itkAffineTransform.h>
#include <itkTransformFileWriter.h>
#include "vtkDoubleArray.h"
#include <vtkPointData.h>
#include <vtkSmartPointer.h>
#include <vtkCenterOfMass.h>

typedef double Point3[3];
using namespace std;

int main( int argc, char * argv[] )
{
	int spharmDegree = 12;  // In our experience, the SPHARM maximal degree is should be chosen between 5 to 50. If the degree is chosen too high the reconstructed SPHARM surface will often show signs of voxelization. (default: 12)
    std::string inParaFile; // input para mesh dataset
    std::string inSurfFile; // input surface mesh dataset
    std::string outSurfFile; // input surface mesh dataset
    std::string bullseyeFile; // output regional curvature in bull's-eye format
	std::string outbase = "./";	// output Directory and base filename
	int x1Position = 0, y1Position = 0, x2Position = 0, y2Position = 0, x3Position = 0, y3Position = 0;
	bool refBasalPoint1=false, refBasalPoint3=false;
	float inSliceSpacing = 0;

	  if( argc != 3 )
	    {
	    std::cerr << "usage: " << argv[0] << " input_filename output_filename" << std::endl;
	    // std::cerr << "  : " << std::endl;
	    exit(1);
	    }

        int subdivLevel = 20; // Subdivision level for linear ikosahedron subdivision (default: 20)
        bool flipTemplateFileOn = false; // Use the Fliptemplate for normalization (default: 0)
	    std::string flipTemplateFile;   // not necessary when flipTemplateFileOn being set to false
	    bool regTemplateFileOn = false; // Use the registration template (default: 0)
	    std::string regTemplateFile;   // not necessary when regTemplateFileOn being set to false
	    bool writePara = false; //The spherical icosahedron subdivision, as well as local phi and theta attribute files for the quality control visualization with KWMeshVisu will be written out. (default: 0)
	    bool NoParaAlignFlag = false; // Do not do any aligment/correspondence of the parametrization using the first order ellipsoid (default: 0)
	    bool debug = false; // Verbose mode (default: 0)
	    bool regParaTemplateFileOn = false; // Use the registration template (default: 0)
	    std::string regParaTemplateFile;   // not necessary when regParaTemplateFileOn being set to false
	    bool regParaPointsOn = false; // Use regParaPoints File (default: 0);
	    std::string regParaPoints;   // not necessary when regParaPointsOn being set to false
	    bool procrustesTransformationOutputOn = false; // Enables the output of the computed Procrustes affine transformation matrix. This can be used to check if improper flips occured. (default: 0)
	    std::string procrustesTransformationFile;     // not necessary when procrustesTransformationOutputOn being set to false
	    int finalFlipIndex = 0; // "This is rarely necessary,only in case the general first order ellipsoid heuristic fails...        Flipping of the parametrization ALONG a given axis will lead to a flip of the coordinates ALONG the same axis in the ellipse aligned mesh.       1 = flip along axes of x & y,     2 = flip along y & z,       3 = flip along x & z       4 = flip along x,        5 = flip along y,        6 = flip along x & y & z,       7 = flip along z  where y is the smallest, x is the second smallest and z is the long axis of the ellipsoid (default: 0)";
	    bool meshWithCurvature = false;


#define dimension 3

  typedef itk::DefaultDynamicMeshTraits<float, dimension, dimension, double> MeshTrait;
  typedef itk::Mesh<float, dimension, MeshTrait>                             MeshType;

  /** Hold on to the type information specified by the template parameters. */
  typedef  MeshType::Pointer    MeshPointer;
  typedef  MeshTrait::PointType MeshPointType;
  typedef  MeshTrait::PixelType MeshPixelType;

  /** Some convenient typedefs. */
  typedef  MeshType::Pointer                MeshPointer;
  typedef  MeshType::CellTraits             CellTraits;
  typedef  MeshType::PointsContainerPointer PointsContainerPointer;
  typedef  MeshType::PointsContainer        PointsContainer;
  typedef  MeshType::CellsContainerPointer  CellsContainerPointer;
  typedef  MeshType::CellsContainer         CellsContainer;
  typedef  MeshType::PointType              PointType;
  typedef  MeshType::CellType               CellType;
  typedef  itk::TriangleCell<CellType>      TriangleType;

  typedef itk::MeshSpatialObject<MeshType>                     MeshSpatialObjectType;
  typedef itk::SpatialObjectWriter<3, float, MeshTrait>        MeshWriterType;
  typedef neurolib::ParametricMeshToSPHARMSpatialObjectFilter  SPHARMFilterType;
  typedef itk::Mesh3DProcrustesAlignFilter<MeshType, MeshType> ProcrustesFilterType;
  typedef itk::AffineTransform<double, dimension>              TransformType;

  std::string  outFileName("dummy");
  const char * base_string;
  base_string = outbase.c_str();


  vtkPolyDataReader *VTKreader = vtkPolyDataReader::New();
  vtkPolyData *      convSurfaceMesh = vtkPolyData::New();

  MeshType::Pointer surfaceMesh = MeshType::New();
  MeshType::Pointer paraMesh = MeshType::New();
  MeshPointer       regParaTemplateMesh;
  MeshPointer       newParaMesh = MeshType::New();
//   typedef itk::SpatialObjectReader<3,float,MeshTrait> ReaderType;
//   ReaderType::Pointer readerSH = ReaderType::New();
  neurolib::SphericalHarmonicSpatialObject::Pointer flipTemplateSO = neurolib::SphericalHarmonicSpatialObject::New();

  try
    {
    if( regParaTemplateFileOn )
      {
      try
        {
        if( debug )
          {
          std::cout << "Reading in regTemplateMesh" << std::endl;
          }

        VTKreader->SetFileName(regParaTemplateFile.c_str() );

        VTKreader->Update();

        vtkPolyData *convTemplateMesh = vtkPolyData::New();
        convTemplateMesh = VTKreader->GetOutput();

        vtkPolyDataToitkMesh* VTKITKConverter4 = new vtkPolyDataToitkMesh;
        VTKITKConverter4->SetInput(convTemplateMesh);
        regParaTemplateMesh = VTKITKConverter4->GetOutput();
        regParaTemplateMesh->Update();

        delete (VTKITKConverter4);

        // 2. use ICP VTK filter -> T
        vtkIterativeClosestPointTransform * ICPTransform = vtkIterativeClosestPointTransform::New();
        vtkLandmarkTransform *              LT = ICPTransform->GetLandmarkTransform();
        LT->SetModeToRigidBody();
        ICPTransform->SetSource(convTemplateMesh);
        ICPTransform->SetTarget(convSurfaceMesh);
        ICPTransform->SetMaximumNumberOfIterations(50);

        ICPTransform->Update();

        }
      catch( itk::ExceptionObject ex )
        {
        std::cout << ex.GetDescription() << "while doing regParaTemplateFile" << std::endl;
        return 1;
        }
      }

    // Main body here
    SPHARMFilterType::Pointer spharmFilter = SPHARMFilterType::New();
    spharmFilter->SetInputSurfaceMesh(surfaceMesh);

    spharmFilter->SetInputParametrizationMesh(paraMesh);// Not necessary

    if( regParaTemplateFileOn )
      {
      spharmFilter->SetInputParametrizationMesh(newParaMesh);
      }
    else
      {
      spharmFilter->SetInputParametrizationMesh(paraMesh);
      }

    spharmFilter->SetDegree(spharmDegree);
    if( NoParaAlignFlag || regParaTemplateFileOn )
      {
      spharmFilter->SetParaEllipseAlignment(false);
      }
    else
      {
      spharmFilter->SetParaEllipseAlignment(true);
      }
    if( flipTemplateFileOn && !regParaTemplateFileOn )
      {
      spharmFilter->SetFlipTemplate(flipTemplateSO);
      }
    if( finalFlipIndex > 0 && finalFlipIndex <= 7 )
      {
      spharmFilter->SetParametrizationFlipIndex(finalFlipIndex);
      }

    clock_t begin, end, start, finish;// by feilang
    time(&begin);
    start = clock();
    spharmFilter->GenerateData(); // **Coefficient Computation
    // spharmFilter->Update();
    time(&end);
    finish = clock();
    std::cout << "GenerateData--Coefficient Computation:" << difftime(end,begin) << std::endl;
    std::cout << "GenerateData--Coefficient Computation:" << finish-start << "--CLOCKS_PER_SEC:" <<CLOCKS_PER_SEC << std::endl;

    if( debug )
      {
      std::cout << "saving par aligned data" << std::endl;
      }

    neurolib::SphericalHarmonicSpatialObject::Pointer      spharmSO = spharmFilter->GetOutput();
    neurolib::SphericalHarmonicSpatialObject::CoefListType coeflist;
    spharmSO->GetCoefs(coeflist); // coefficients stored in coeflist

    // save associated surface
    neurolib::SphericalHarmonicMeshSource::Pointer meshsrc = neurolib::SphericalHarmonicMeshSource::New();
    meshsrc->SetCoefs(coeflist);
    meshsrc->SetLevel(subdivLevel);
    meshsrc->Update();

    MeshType* meshSH;
    meshSH = meshsrc->GetOutput();

    itkMeshTovtkPolyData * ITKVTKConverter = new itkMeshTovtkPolyData;
    ITKVTKConverter->SetInput( meshSH );

    vtkPolyDataWriter *vtkwriter;
    vtkwriter = vtkPolyDataWriter::New();

    meshWithCurvature = true;

    bool bullEye=true;
    if(bullEye)
    {
    	// Compute the center of mass
    	vtkSmartPointer < vtkCenterOfMass > centerOfMassFilter = vtkSmartPointer< vtkCenterOfMass > ::New();
    	centerOfMassFilter->SetInput(ITKVTKConverter->GetOutput() );
    	centerOfMassFilter->SetUseScalarsAsWeights(false);
    	centerOfMassFilter->Update();
    	double center1[3];
    	centerOfMassFilter->GetCenter(center1);
    	std::cout << "Mass center is: " << center1[0] << " " << center1[1]<< " " << center1[2] << std::endl;

    	float split1_2 = M_PI*5/3;
    	if(refBasalPoint1||refBasalPoint3)
    	{
    		split1_2 = atan2(-x1Position*inSliceSpacing-center1[0], -y1Position*inSliceSpacing-center1[1]); // in the interval [-pi,+pi] radians
    		split1_2 = split1_2+M_PI; //in the interval [0,+2pi] radians
    	}
    	std::cout << "split1_2 is: " << split1_2*180/M_PI << " degrees." << std::endl;
    	float split2_3 = split1_2 + M_PI/3;
    	split2_3 = (split2_3<=2*M_PI) ? split2_3 : (split2_3-2*M_PI);
    	float split3_4 = split2_3 + M_PI/3;
    	split3_4 = (split3_4<=2*M_PI) ? split3_4 : (split3_4-2*M_PI);
    	float split4_5 = split3_4 + M_PI/3;
    	split4_5 = (split4_5<=2*M_PI) ? split4_5 : (split4_5-2*M_PI);
    	float split5_6 = split4_5 + M_PI/3;
    	split5_6 = (split5_6<=2*M_PI) ? split5_6 : (split5_6-2*M_PI);
    	float split6_1 = split5_6 + M_PI/3;
    	split6_1 = (split6_1<=2*M_PI) ? split6_1 : (split6_1-2*M_PI);

    	float split7_8 = split1_2;
    	if(refBasalPoint3)
    	{
        	split7_8 = atan2(-x2Position*inSliceSpacing-center1[0], -y2Position*inSliceSpacing-center1[1]); // in the interval [-pi,+pi] radians
        	split7_8 = split7_8+M_PI; //in the interval [0,+2pi] radians
    	}
    	std::cout << "split7_8 is: " << split7_8*180/M_PI << " degrees." << std::endl;
    	float split8_9 = split7_8 + M_PI/3;
    	split8_9 = (split8_9<=2*M_PI) ? split8_9 : (split8_9-2*M_PI);
    	float split9_10 = split8_9 + M_PI/3;
    	split9_10 = (split9_10<=2*M_PI) ? split9_10 : (split9_10-2*M_PI);
    	float split10_11 = split9_10 + M_PI/3;
    	split10_11 = (split10_11<=2*M_PI) ? split10_11 : (split10_11-2*M_PI);
    	float split11_12 = split10_11 + M_PI/3;
    	split11_12 = (split11_12<=2*M_PI) ? split11_12 : (split11_12-2*M_PI);
    	float split12_7 = split11_12 + M_PI/3;
    	split12_7 = (split12_7<=2*M_PI) ? split12_7 : (split12_7-2*M_PI);

    	float split13_14 = split1_2 + M_PI/12;
    	if(refBasalPoint3)
    	{
        	split13_14 = atan2(-x3Position*inSliceSpacing-center1[0], -y3Position*inSliceSpacing-center1[1]); // in the interval [-pi,+pi] radians
        	split13_14 = split13_14+M_PI; //in the interval [0,+2pi] radians
    	}
    	std::cout << "split13_14 is: " << split13_14*180/M_PI << " degrees." << std::endl;
    	float split14_15 = split13_14 + M_PI_2;
    	split14_15 = (split14_15<=2*M_PI) ? split14_15 : (split14_15-2*M_PI);
    	float split15_16 = split14_15 + M_PI_2;
    	split15_16 = (split15_16<=2*M_PI) ? split15_16 : (split15_16-2*M_PI);
    	float split16_13 = split15_16 + M_PI_2;
    	split16_13 = (split16_13<=2*M_PI) ? split16_13 : (split16_13-2*M_PI);


        MeshType *             paraMesh = meshsrc->GetOutputParaMesh();
        PointsContainerPointer paraPoints = paraMesh->GetPoints();
        int nvert = paraPoints->Size();
    	double * phi_theta = meshsrc->GetPrecomputedIcosList(); // nvert*2
    	double * phi = new double[nvert];
    	double * theta = new double[nvert];
    	double * curvedness = meshsrc->GetCurvedness();         // nvert*1

    	double curvednessBullsEye[17];
    	for(int i=0; i<17; ++i)
    		curvednessBullsEye[i]=0;

    	int nvertBullEye[17]; // number of vertex in each bull eye
    	for(int i=0; i<17; ++i)
    		nvertBullEye[i]=0;
// HYL
    	//outFileName.erase();
    	//outFileName.append("MyCurvature.txt");
    	//std::ofstream myefile(outFileName.c_str(), std::ios::out);
    	//myefile.precision(10);
    	//for (int i=0; i<nvert; i++)
    	//{
    	//	    myefile << "Curvedness["<<i<<"]: " << curvedness[i] << std::endl;
    	//}
    	//myefile.close();

    	//double thickness[nvert];
    	//std::ofstream myefile;
    	//myefile.open ("Thickness.txt", std::ios::in);

    	FILE *fd;
    	if(0==(fd=fopen(argv[1], "r"))) {


    	printf("error in reading file\n");
    	}
    	char param_name[200], param_str[200];
    	float param[nvert];
    	double param_db[nvert];
    	for (int i=0; i<nvert; i++)
    	{
    		fscanf(fd,"%s %f", param_name, &param[i]);
    		printf("%s=%f\n",param_name,param[i]);
    		param_db[i]=param[i];
    	}
    	fclose(fd);

    	//for (int i=0; i<nvert; i++)
    	//{
    	//	myefile >> thickness[i];
       	//	cout << "Thickness["<<i<<"]: " << thickness[i] << std::endl;
    	//}
    	curvedness = param_db;

// HYL //


    	for( int i = 0; i < nvert; i++ )
    	{
    		phi[i] = phi_theta[2 * i]; // 0~2pi
    		theta[i] = phi_theta[2 * i + 1]; // 0~pi


    	if( (theta[i]>=0) && (theta[i]<M_PI/3) && ( ((split6_1<split1_2)&&(phi[i]>=split6_1)&&(phi[i]<split1_2)) || ((split6_1>split1_2)&&((phi[i]>=split6_1)||(phi[i]<split1_2))) ) ) // 0~60;0~60;
    	{
    		curvednessBullsEye[1] += curvedness[i];
    		nvertBullEye[1]++;
    	}
    	else if( (theta[i]>=0) && (theta[i]<M_PI/3) && ( ((split1_2<split2_3)&&(phi[i]>=split1_2)&&(phi[i]<split2_3)) || ((split1_2>split2_3)&&((phi[i]>=split1_2)||(phi[i]<split2_3))) ) ) // 0~60;60~120;
    	{
    		curvednessBullsEye[2] += curvedness[i];
    		nvertBullEye[2]++;
    	}
    	else if( (theta[i]>=0) && (theta[i]<M_PI/3) && ( ((split2_3<split3_4)&&(phi[i]>=split2_3)&&(phi[i]<split3_4)) || ((split2_3>split3_4)&&((phi[i]>=split2_3)||(phi[i]<split3_4))) ) ) // 0~60;120~180;
    	{
    		curvednessBullsEye[3] += curvedness[i];
    		nvertBullEye[3]++;
    	}
    	else if( (theta[i]>=0) && (theta[i]<M_PI/3) && ( ((split3_4<split4_5)&&(phi[i]>=split3_4)&&(phi[i]<split4_5)) || ((split3_4>split4_5)&&((phi[i]>=split3_4)||(phi[i]<split4_5))) ) ) // 0~60;180~240;
    	{
    		curvednessBullsEye[4] += curvedness[i];
    		nvertBullEye[4]++;
    	}
    	else if( (theta[i]>=0) && (theta[i]<M_PI/3) && ( ((split4_5<split5_6)&&(phi[i]>=split4_5)&&(phi[i]<split5_6)) || ((split4_5>split5_6)&&((phi[i]>=split4_5)||(phi[i]<split5_6))) ) ) // 0~60;240~300;
    	{
    		curvednessBullsEye[5] += curvedness[i];
    		nvertBullEye[5]++;
    	}
    	else if( (theta[i]>=0) && (theta[i]<M_PI/3) && ( ((split5_6<split6_1)&&(phi[i]>=split5_6)&&(phi[i]<split6_1)) || ((split5_6>split6_1)&&((phi[i]>=split5_6)||(phi[i]<split6_1))) ) ) // 0~60;300~360;
    	{
    		curvednessBullsEye[6] += curvedness[i];
    		nvertBullEye[6]++;
    	}
    	else if( (theta[i]>=M_PI/3) && (theta[i]<M_PI*2/3) && ( ((split12_7<split7_8)&&(phi[i]>=split12_7)&&(phi[i]<split7_8)) || ((split12_7>split7_8)&&((phi[i]>=split12_7)||(phi[i]<split7_8))) ) ) // 60~120;0~60;
    	{
    		curvednessBullsEye[7] += curvedness[i];
    		nvertBullEye[7]++;
    	}
    	else if( (theta[i]>=M_PI/3) && (theta[i]<M_PI*2/3) && ( ((split7_8<split8_9)&&(phi[i]>=split7_8)&&(phi[i]<split8_9)) || ((split7_8>split8_9)&&((phi[i]>=split7_8)||(phi[i]<split8_9))) ) ) // 60~120;60~120;
    	{
    		curvednessBullsEye[8] += curvedness[i];
    		nvertBullEye[8]++;
    	}
    	else if( (theta[i]>=M_PI/3) && (theta[i]<M_PI*2/3) && ( ((split8_9<split9_10)&&(phi[i]>=split8_9)&&(phi[i]<split9_10)) || ((split8_9>split9_10)&&((phi[i]>=split8_9)||(phi[i]<split9_10))) ) ) // 60~120;120~180;
    	{
    		curvednessBullsEye[9] += curvedness[i];
    		nvertBullEye[9]++;
    	}
    	else if( (theta[i]>=M_PI/3) && (theta[i]<M_PI*2/3) && ( ((split9_10<split10_11)&&(phi[i]>=split9_10)&&(phi[i]<split10_11)) || ((split9_10>split10_11)&&((phi[i]>=split9_10)||(phi[i]<split10_11))) ) ) // 60~120;180~240;
    	{
    		curvednessBullsEye[10] += curvedness[i];
    		nvertBullEye[10]++;
    	}
    	else if( (theta[i]>=M_PI/3) && (theta[i]<M_PI*2/3) && ( ((split10_11<split11_12)&&(phi[i]>=split10_11)&&(phi[i]<split11_12)) || ((split10_11>split11_12)&&((phi[i]>=split10_11)||(phi[i]<split11_12))) ) ) // 60~120;240~300;
    	{
    		curvednessBullsEye[11] += curvedness[i];
    		nvertBullEye[11]++;
    	}
    	else if( (theta[i]>=M_PI/3) && (theta[i]<M_PI*2/3) && ( ((split11_12<split12_7)&&(phi[i]>=split11_12)&&(phi[i]<split12_7)) || ((split11_12>split12_7)&&((phi[i]>=split11_12)||(phi[i]<split12_7))) ) ) // 60~120;300~360;
    	{
    		curvednessBullsEye[12] += curvedness[i];
    		nvertBullEye[12]++;
    	}
    	else if( (theta[i]>=M_PI*2/3) && (theta[i]<=M_PI) && ( ((split16_13<split13_14)&&(phi[i]>=split16_13)&&(phi[i]<split13_14)) || ((split16_13>split13_14)&&((phi[i]>=split16_13)||(phi[i]<split13_14))) ) ) // 120~180;-15(345)~75;
    	{
    		curvednessBullsEye[13] += curvedness[i];
    		nvertBullEye[13]++;
    	}
    	else if( (theta[i]>=M_PI*2/3) && (theta[i]<=M_PI) && ( ((split13_14<split14_15)&&(phi[i]>=split13_14)&&(phi[i]<split14_15)) || ((split13_14>split14_15)&&((phi[i]>=split13_14)||(phi[i]<split14_15))) ) ) // 120~180;75~165;
    	{
    		curvednessBullsEye[14] += curvedness[i];
    		nvertBullEye[14]++;
    	}
    	else if( (theta[i]>=M_PI*2/3) && (theta[i]<=M_PI) && ( ((split14_15<split15_16)&&(phi[i]>=split14_15)&&(phi[i]<split15_16)) || ((split14_15>split15_16)&&((phi[i]>=split14_15)||(phi[i]<split15_16))) ) ) // 120~180;165~255;
    	{
    		curvednessBullsEye[15] += curvedness[i];
    		nvertBullEye[15]++;
    	}
    	else if( (theta[i]>=M_PI*2/3) && (theta[i]<=M_PI) && ( ((split15_16<split16_13)&&(phi[i]>=split15_16)&&(phi[i]<split16_13)) || ((split15_16>split16_13)&&((phi[i]>=split15_16)||(phi[i]<split16_13))) ) ) // 120~180;255~345;
    	{
    		curvednessBullsEye[16] += curvedness[i];
    		nvertBullEye[16]++;
    	}
    	else
    	{
    		std::cout << "theta[i] = "<< theta[i] << ";  phi[i] = "<< phi[i] << std::endl;
    		std::cout << "split16_13 = "<< split16_13 << ";  split13_14 = "<< split13_14 << std::endl;
    		std::cout << "split14_15 = "<< split14_15 << ";  split15_16 = "<< split15_16 << std::endl;
    	}
    	}
    	std::cout << std::endl;

        outFileName.erase();
        outFileName.append(bullseyeFile);
 //       std::ofstream efile("Thickness_Bullseye.txt", std::ios::out);
        std::ofstream efile(argv[2], std::ios::out);
    	   efile.precision(10);
    	   efile << "NUMBER_OF_SEGMENTS = " << 16 << std::endl;

    	for (int i=1; i<=16; ++i)
    	{
    		curvednessBullsEye[i] = curvednessBullsEye[i]/nvertBullEye[i];
    		std::cout << "nvertBullEye["<<i<<"]: " <<nvertBullEye[i]<< ";  ThicknessBullsEye["<<i<<"]: "<<curvednessBullsEye[i]<< std::endl;
    		efile << curvednessBullsEye[i] << std::endl;
    	}
    	efile.close();
    }

    }

  catch( itk::ExceptionObject e )
    {
    e.Print(std::cout);
    return EXIT_FAILURE;
    }

  return EXIT_SUCCESS;

}
